#include "StaticSolver.h"

#include <iostream>

#include "Utility/Logging/SGLogger.h"
#include "Utility/Algebra/Factory.h"
#include "Utility/Algebra/LinearSolver.h"

#ifdef SG_LICENSE_CHECK
#include "Utility/Security/LicenseCheck.h"
#endif

#include "DBManager/IService/Input/Subcase/ISubcaseService.h"
#include "DBManager/IService/Input/Element/IElementService.h"
#include "DBManager/IService/Input/Load/ILoadCombinationService.h"
#include "DBManager/IService/Result/Nodal/IDISPLACEMENTService.h"

#include "SGFem/Pre/Assembly/Assembly.h"
#include "SGFem/Pre/Assembly/DataAdapter.h"

#include "SGFem/ElementCalculator/Factory/ElementCalculatorFactory.h"
#include "SGFem/Task/Kernel/MModel.h"

using SG::DataStructure::Common::Status;
using SG::DataStructure::Common::Id;
using SG::DataStructure::Common::IdArray;
using SG::DataStructure::FEM::Project;
using SG::DataStructure::FEM::Subcase;
using SG::DataStructure::FEM::LoadCombination;
using SG::DataStructure::FEM::SOLType;

using namespace SG::DBManager;

using SG::FEMSolver::EigenVector;
using SG::FEMSolver::AssemblyEigen;
using SG::FEMSolver::Element::ElementCalculatorFactory;

class StaticSolver::Impl
{
public:
    Impl();
    Impl (const Impl&) = delete;
    Impl& operator= (const Impl&) = delete;
    ~Impl ();
    std::shared_ptr<SG::Algebra::LinearSolver> s_solver;
};

StaticSolver::Impl::Impl ()
{
    // 初始化稀疏矩阵求解器
    // std::string solver_type = "sgEigenLLT";
    // std::string solver_type = "sgSuperLU";
    // std::string solver_type = "sgPetscKsp";
    std::string solver_type = "TPardiso<Real_t>";
    auto        solver      = reinterpret_cast<SG::Algebra::LinearSolver*> (SG::Base::Factory::GetInstance ()->Produce (solver_type));
    if (solver)
    {
        s_solver.reset (solver);
    }
    assert (solver && "Make sure that the required solver exists!");
}

StaticSolver::Impl::~Impl ()
{
    // perf: force to release the resources allocated for the linear solver.
    s_solver.reset ((SG::Algebra::LinearSolver* ) (nullptr));
}

StaticSolver::StaticSolver() : m_impl(new Impl())
{
}

StaticSolver::~StaticSolver()
{
   delete m_impl;
}

void StaticSolver::Solve (const Kernel::MModel& model, const Id loadId, _OUT std::vector<SG::DataStructure::Common::Real>& displacement)
{
    // 多工况 多工况处理可能在单元计算之后处理
    // 创建分析类型，分析
    auto& assembly = AssemblyEigen::GetInstance();

    EigenVector f;
    model.ComputeLoadVector(loadId, f);

    SGTrace << "back substitution start";
    EigenVector displacement0;
    m_impl -> s_solver->solve (f, _OUT displacement0);
    if (SG::Algebra::Info_t::Success != m_impl -> s_solver -> info ())
    {
        std::cerr << "solve failed !!!" << std::endl;
        // return;
    }
    SGTrace << "back substitution fini";

    assembly.RecoveryDisplacement (displacement0, _OUT displacement);
    std::cout << "== solve Fini"<< std::endl;
    // 积分点应力应变计算， 外推插值

    // 应力磨平

    std::cout << "static Solve fini  \n";
}

void StaticSolver::Solve (const SG::DBManager::DBServiceFactorySP dbServiceFactory, const SG::DataStructure::FEM::Project& project, const SG::DataStructure::Common::IdArray& ids)
{
    // 检查授权
    #ifdef SG_LICENSE_CHECK
        SG::Utility::SGCheckLicense();
    #endif
    
    Kernel::MModel model (SOLType::SOL_SESTATIC, dbServiceFactory);
    model.Initialize (ids);
    model.CalcMatrixEquation (ids[0]);

    SG::Algebra::SpMat_t K;
    model.ComputeStiffnessMatrix (_OUT K);

    SGTrace << "decomp start";
    m_impl -> s_solver -> compute(K);
    if (SG::Algebra::Info_t::Success != m_impl -> s_solver -> info ())
    {
        std::cerr << "Decomposition failed !!!" << "fail type " << m_impl -> s_solver -> info () << std::endl;
        return;
    }
    std::cout << "== Matrix Decomposition" << std::endl;
    SGTrace << "decomp fnit";

    // 依次计算每个工况的位移
    const auto& dispService    = dbServiceFactory->get<IDISPLACEMENTService> ();
    const auto& subcaseService  = dbServiceFactory -> get<ISubcaseService>();
    for(auto subcaseId : ids)
    {   
        // 获取LOAD的荷载集ID
        auto para = subcaseService -> FindSubCase (subcaseId);
        Subcase subcase(para);
        std::cout << "== Get RHS " << std::endl;
        // 计算位移
        std::vector<SG::DataStructure::Common::Real> displacement;
        StaticSolver::Solve (model, subcase.GetLoadId (), _OUT displacement);

        // 保存静力工况位移
        std::vector<SG::DataStructure::Post::DISPLACEMENTData> displacementData;
        SG::FEMSolver::DataAdapter::PackGridField (displacement, _OUT displacementData);
        dispService->Add ({SG::DataStructure::FEM::SOLType::SOL_SESTATIC, subcaseId}, displacementData);

        // 积分点应力应变计算， 外推插值
        std::cout << "===== compute Integ stress \n";
        auto elementTypes = dbServiceFactory->get<IElementService> ()->GetAllElementTypes ();
        for (const auto& type : elementTypes)
        {
            auto eleCalculator = ElementCalculatorFactory::Create (type);
            if (NULL == eleCalculator)
            {
                std::cout << "can 't find element calc "<< std::endl;
            }
            
            eleCalculator -> Initialize (type);
            eleCalculator -> DataRecovery (dbServiceFactory, subcaseId, 0.0, displacement);
        }
    }
}